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ABSTRACT 



Military electronics are frequently operated in excessively confined spaces aboard 
ships and aircraft. This limited space impacts significantly on the space available 
for cooling equipment. The optimal solution is the development of one universal, 
modular rack for shipboard or aviation use. With a modular design, upgrades to 
equipment could also be accompanied by an upgrade to the cooling rack itself with 
very little additional cost or difficulty. A water cooled avionics rack can provide suf- 
ficient cooling for any piece or combination of avionics equipment if enough water 
flow paths are used, the water is at the appropriated temperature and the water 
is properly distributed within the passages. To determine if the cooling medium, 
water, is sufficiently distributed within a modular cooling rack, an analysis of the 
flow and pressure distribution of the coolant is required. This thesis presents a 
computer code that has been developed as an initial step in the total design of a 
modular cooling rack for avionics equipment. In itself, the code details a specific 
design technique and allows for the determination of whether a proposed configu- 
ration, including source location, characteristics of the cooling water, and the size 
and shape of the proposed flow passages will indeed provided a proper distribution 
of the coolant. 
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I. INTRODUCTION 



A. GENERAL 

This thesis describes the development of a computer code to analyze the flow 
and pressure distribution of water in an avionics cooling rack. This code allows for 
variable rack design and flow sources. The code is written in FORTRAN 77. In 
addition, the code can be run on any IBM or IBM compatible personal computer. 

B. BACKGROUND 

Extremely high temperatures are a primary cause of avionics equipment un- 
reliability. The origin of the thermal problem is in the continuous effort to develop 
lighter and smaller components. These smaller, more densely packed components, 
by virtue of the heat dissipation within a small volume, frequently operate at exces- 
sively high temperatures. These high temperatures result in component derating, 
performance degradation and accelerated failure. [Ref 1] 

Successful thermal management of electronic systems, under development in 
the latter part of the 1990s will require the removal of as much as 500 W from a 
single chip, at heat fluxes between 50 to 100 W/cm^. Volumetric heat release rates 
can also be expected to increase dramatically and are likely to exceed 10 W/cm^. 
Silicon chips, with embedded bipolar circuits, have traditionally been maintained at 
temperatures ranging from 65° C, for commercial computers, to between 110° C and 
125° C for military equipment. [Ref 2] 

Due to necessity, military electronics are frequently operated in confined 
spaces aboard ships and aircraft. This limited space impacts significantly on the 
space available for extensive cooling equipment. Air-cooled avionics systems are 
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currently the most frequently encountered cooling systems used in military appli- 
cations. However, it is common practice with military avionics to upgrade the 
operational capability of a ship or aircraft by upgrading only one or two pieces of 
avionics equipment at a time. With an air cooled system, this is very ineffective. 
By replacing equipment in one of these systems with equipment of a different size 
and shape, the cooling airflow of the original design is disturbed and the resulting 
airflow around the new (as well as the old) equipment is less than optimal. The 
original cooling system was designed for equipment configured in a certain manner 
and little attention is paid to the reoptimization of the cooling system when changes 
are made. Moreover, cost restrictions are also an important factor in the design of 
military cooling systems. 

One solution for the military applications problem is the development of a 
universal rack for shipboard or aviation use. If this rack were modular, it would allow 
for a great amount of flexibility in design and provide significant cost savings. With 
a modular design, upgrades to equipment could also be accompanied by an upgrade 
to the cooling rack itself with very little additional cost or difficulty. By developing 
one universal rack for all military applications, a significant improvement in cooling 
systems for updated designs can also be realized. This could be accompanied by 
additional cost savings through the elimination of unique cooling systems for every 
different avionics suite. 

The rack would be a structure that could accomodate the placement of several 
modular electronic packages. One might imagine a ^shoebox^ configuration with, say, 
16 ‘‘/lo/es” in which 16 electronic packages could be inserted. Each of the packages is 
presumed to contain electronic components mounted in a variety of ways (i.e printed 
circuit boards or on brackets attached to the walls of the package). Such a structure 
would possess at least 64 flow passages exclusive of flow passages that are part of 
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the fluid source or pump. The objective is to assure that each flow passage carries 
enough fluid to absorb the dissipated heat that is somehow conducted to the rack 
structure. Of course, a heat exchanger may be required to transfer the heat to the 
ultimate heat sink which may be the environmental air or sea water. 

Thus, the bcisic assumption is made that a modular, water cooled avionics 
rack can provide sufflcient cooling for any piece or combination of avionics equipment 
if enough water flow paths are used, the water is at the appropriate temperature and 
the water is properly distributed within the passages. To determine if the cooling 
medium, water, is being sufficiently distributed within a modular cooling rack, an 
analysis of the flow and pressure distribution of the water is required. This can be 
accomplished using a computer code utilizing node analysis. Such a computer code, 
allowing for a variable number of branches and junctions, is presented here as a first 
step in the development of a universal cooling rack for military applications. 

C. BASIC THEORY AND APPROACH 

1. Node Analysis 

The analysis of the cooling rack is based on the stipulation that any 
size (diameter and length) of passage may be used in the construction of the rack. 
Variables in the program include the density and viscosity of the water, the number 
and location of flow sources, the number of water paths entering each junction, the 
shape of the flow passages (circular or rectangular) and the length and equivalent 
diameter of each passage. The purpose of this degree of flexibility is to allow for easy 
redesign of a rack in the event that it does not meet the requirements which are the 
proper distribution of cooling water in each section of pipe. By varying either the 
rack configuration or the state of the water via computer input, a rack that provides 
the proper flow distribution can eventually be proposed. 
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The calculations for the flow in the passages employs a matrix oriented 
procedure used in network analysis. The network analysis approach can be tailored 
to flow in paissages by proposing an analogy of the ’’current” to fluid flow and the 
"voltage” to the fluid pressure. The "resistance” is then attributed to the friction in 
each length of pipe. Therefore, each length of pipe will have a resistance associated 
with it, and possibly a pressure source as well, depending on the design of the cooling 
rack. 

2. Laminar Flow 

The computer code is designed to calculate a laminar flow distribution. 
If the flow is in transition or turbulent, there is a significant increase in the amount 
of frictional resistance. For Reynolds numbers (Re) less than or equal to 2100, the 
flow is considered laminar. For Reynolds number between 2100 and 10000, the flow 
is in transition and the flow is turbulent if the Reynolds number exceeds 10000. The 
computer code calculates the Reynolds number for each flow passage and provides 
a warning if the Reynolds number exceeds 2100. 

D. SCOPE 

• Chapter II explains and details the basic code required to analyze the flow 
and pressure distribution of water in a proposed cooling rack design. 

• Chapter III describes the computer code, its essential capabilities and limita- 
tions, associated subroutines, input requirements and final output. 

• Chapter IV presents the results from several case runs that exhibit the flexi- 
bility and capability of the method. 

• Chapter V concludes with future development efforts and the application po- 
tential of this computer code. 
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II. AVIONICS COOLING RACK PROBLEM 

FORMULATION 

A. GENERAL 

The flow distribution of the fluid (in all branches) of the cooling rack is 
determined using a matrix oriented solution technique. The general equations and 
the matrix solution are presented in what follows. 

1. COOLING RACK BASIC STRUCTURE 

The basic structure of the cooling rack is a series of fluid passages fit- 
ted together using standard junction components. An almost unlimited number of 
passages may be fitted together. An example of one possible simple combination of 
passages with one flow source is presented in Figure 2.1. 




This example illustrates how variable sized flow passages may be used. 
The rack is composed of b straight passages called branches and nt junctions called 
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nodes. The branches are numbered because the computer code requires that basic 
information for each flow passage branch be entered separately. The nodes are also 
numbered and those are shown in circles. Note also that there is a pressure source 
located at branch- 1. For convenience, the numbering sequence in Figure 2.1, begins 
with a 1 on the upper left side and proceeds to the right and then down. Although 
this numbering sequence is arbitrary, maintaining a consistent technique facilitates 
later adaptations to the rack and comparison with other rack configurations. 

B. VARIABLES IN DESIGN 

As stated in a previous section, the length, shape and the diameters of the flow 
passages may vary. The code is written for either circular or rectangular passages. 
For rectangular passages an effective diameter is calculated. Although the actual 
diameter or effective diameter may vary within one rack design, it is assumed that 
either all circular or all rectangular passages are used for any one rack configuration. 
All length and diameter size information is entered in the initial part of the program 
as 6 X 1 matrices named ELL (IB) and D(IB), respectively. The density and the 
viscosity of the cooling water are input as the variables RHO and MU. The other 
input and output variables used in the computer code are summarized, for the 
reader’s convenience, in Table 1. 



C. DEVELOPMENT AND LINEARIZATION OF 
FLOW EQUATIONS 

1. Laminar Equations 

The basic equation used to determine the change in pressure or pressure 
loss within a fluid passage is derived from the D’Arcy equation 



hj 



fLV^ 

2gdc 



(2.1) 
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TABLE 2.1: INPUT AND OUTPUT VARIABLES (FORTRAN) 



Variable 


Explanation 


A1 


Cross sectional height 


B 


Number of branches 


B1 


Cross sectional width 


D 


Branch diameter vector 


ELL 


Branch length vector 


IB 


Counter for branches, IB = 


MU 


Viscosity 


N 


Number of nodes 


NF 


Node “from” array 


NT 


Node “to” array 


PS 


Pressure source vector 


RHO 


Density 



In this equation, hf \s the pressure loss due to friction and is the 
equivalent diameter. The D’Arcy equation is valid for steady flow within passages 
running full of liquid. In eq (2.1), / is a friction factor. If the flow is laminar, the 
friction factor, /, as shown in Figure 2.2, can be represented by the equation 




( 2 . 2 ) 



The Reynolds number, by definition is 



Therefore, /, is seen to have an alternate form 

^ 64/i 
^ pVd, 

Substitution of eq (2.4) into eq (2.1) yields 



h 



/ = 



S2LVn 

pgdl 



(2.3) 



(2.4) 



(2.5) 
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Friellon FoctOf, 




Figure 2.2; Moody’s chart for the friction factor. 
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By definition, the flow rate, Q, is 



Q = VA (2.6) 

The area. A, of a circular flow passage is, of course, given by 

A = (2.7a) 

but for rectangular passages of side dimensions a and 6, it is 

A = ab (2.76) 

and for square passages where a = 6, it is 

A = a'^ (2.7c) 

For circular passages, dg, is simply the passage diameter, d, and for rectangular 
passages, d^ is defined as 

CA 



— p 

where A is the passage flow area and P, for channels flowing full, is the pcissage 
wetted perimeter. In order to make the equivalent diameter for a circular passage 
equal to the actual diameter, d,C = 4 so that for a rectangular passage with side 
dimensions, a and 6, the equivalent diameter becomes 

2ab 



de — , 

(2 “h 0 



(2.8a) 



In the event that the passage is square (a = 6), the equivalent diameter becomes 



d« a 



( 2 . 86 ) 



Substitution of eqs (2.7) and (2.8) into equation (2.5) yields for the 
circular passage 



hf = 



\28LQfi 

gwd'^p 



(2.9a) 
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for the rectangular passage 



h f = 



SLfiQja + b)^ 
pg{ab)^ 



and for the square passage 



, Z2LpQ 

b/ = — 

pga* 



(2.96) 



(2.9c) 



To simplify eq (2.9a) further, let 

hj = P 



and 



R = 



mip 



This yields the equation 



which is of the same form as 



gird^p 

P=RQ 
V = RI 



( 2 . 10 ) 



Similar adjustments can be made to eqs (2.9b) and (2.9c). When these adjustments 
are made, it is seen that for a rectangular passage, eq (2.9b) gives 

8Lp{a + 6 )^ 



^ P5(a6)3 

and that for a square passage eq (2.9c) provides 

32Lp 



( 2 . 11 ) 



R = 



pga 



Additional losses (other than passage friction already discussed) occur 
in flow passage systems and cannot be disregarded without appreciable error. Com- 
pensation for entrance and exit losses must be considered in the cooling rack system 
by adding an equivalent length of straight pipe. The equivalent length used for 
square edged entries is 20 diameters (or equiv'alent diameters) and for exits, the 
equivalent length is 40 diameters (or equivalent diameters) [Ref. 2]. Therefore, the 
additional length of 60 passage diameters accounts for both exit and entrance losses. 
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D. NUMERICAL SOLUTION SCHEME 

We first consider a general branch (the branch) as shown in Figure 2.3. 
The branch can contain a pump with pressure head, p,fc, with an associated section 
of pipe of resistance, Rk. The branch carries a flow, g*, and possesses a total head 
loss, pk- Notice that the orientation of the head loss is in the direction of the flow 
(from to 




Also observe the orientation of the flow, 9 *, and notice that by the analogy 
to Ohm’s law and because of compatibility (the sum of the losses must match the 
head available) 

Pk = Rk9k+ Psk 

This may be written for all k branches in matrix form 

P = RQ + P, 

Consider a pump and a length of discharge piping. The pump can be repre- 
sented as a flow source or as a pressure source as indicated in Figure 2.4. 
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Figure 2.4: Alternative source arrangements for the development of the 
Flow source Pressure source transformation. 



In these figures, the R's represent the resistance to flow of the discharge 
piping and, in the case of the pressure source, the + and — indicate the discharge 
(-f ) and suction ( — ). In Figure 2.4a, continuity dictates that, 



9 = 




or 



P — d” 



In Figure 2.4b, compatability (the pressure drops around the simple fluid loop must 
match) 

p = Rpp, 



If there is to be an equivalence between the two (Figures 2.4a and 2.4b) then 



R = Rp = Rq 



and 



Ps = Rq, 
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or indeed 




Thus, when we represent a pump as a pressure source with a resistance, we can 
immediately transform this to an equivalent flow source. The computer code is 
written for the source input to be a flow source. 

Note also that continuity dictates that for the branch 

Qk ~ ^kPk d" 

And for all k branches this can be written in matrix form as 



Q = YP + Qs 

Next consider a rack containing b branches and n, nodes. The n\^ node is the 
datum node (the node at the suction end of the pump), and we may set down an 
Tit X b augmented node-branch incidence matrix A“ with elements 



Oy = < 



' -f-1 if branch j leaves node i 
— 1 if branch j enters node i 
0 if branch j is not incident 
upon node i 



( 2 . 10 ) 



For example in the network displayed in Figure 2.5 with its oriented graph shown 
in Figure 2.6, there are = 5 nodes and 6 = 6 branches. For this network, the 
augmented node branch incidence matrix A“ is 



A“ = 



1 

0 

0 

0 

-1 



1 

-1 

0 

0 

0 



0 0 

1 1 

-1 0 

0 -1 

0 0 



0 0 

0 0 

1 0 

-1 1 

0 -1 



Note that for each column, which represents a single branch, there should 
only be two non-zero entries, a 1 and a —1. These entries correspond to the nodes 
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Figure 2.5: Network for example. 




Figure 2.0: Oriented graph for network in Figure 2.5. 
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at the extremities of the branch. A branch that leaves a node shows a +1 while a 
branch that leaves a node shows a —1. 



Next, we form a node-branch incidence matrix which is n x 6. This is done 
by deleting the row in A“ that corresponds to the datum node. In our example, 
row-5 is the datum node and thus 



A = 



1 1 0 

0 -1 1 

0 0-1 

0 0 0 



0 0 0 

1 0 0 

0 1 0 

-1 -1 1 



The vector Q is defined as a 6 x 1 vector representing the flow in all k branches 
and its elements are designated as qk- The product of A with Q, AQ, represents a 
series of equations, one equation for each node. These equations represent the sum 
of all flow in and out of the respective node. Continuity demands that product, AQ 
be null 

AQ = 0 

and this simple matrix equation is a statement of continuity for the whole system. 
To illustrate 



AQ = 



1 1 0 0 0 0 

0-1 1 1 00 

0 0-1 0 10 

0 0 0 -1 -1 1 



91 

92 




{91 + 92 ) 


93 




i~92 + 93 + 94) 


94 




(“93 + 9s) 


95 




. (“94 “ 9s + 9s) . 


. 96 . 









■ 0 ■ 




0 




0 




1 

0 

1 



The pressures or heads at the nodes are represented by an n x 1 vector 
designated by H and the branch pressure drops are represented by an n x 1 vector 
designated by P. Reference to Figure 2.6 shows that 



Pi = hi 
?2 = hi — /l2 



15 



Pz — h2 — hz 

P4 = hz — h.4 

Ps = hz — /14 

PS ^4 



The matrices P and H can be related by a matrix C 



P = CH = 



and it can be observed that 



1 



0 



1 -1 



0 

0 

0 

0 



0 

0 



0 
0 
0 

0 -1 

1 -1 



1 -1 

1 
0 

0 0 1 



hi 

hz 

hs 

h4 



C = A 



Refer now to the general branch in Figure 2.3 and the vector representation 
for the flow in all k branches. The 6x1 vector, Q which represents the branch flow 
rates has been shown to be 

Q = YP + Qs (2.11) 



Again, we carefully note the orientation of the flow source. Here Y is an n x 
n source admittance matrix. Premultiply eq (2.11) by A to obtain 



AQ = AYP + AQs 



But, P = A^H and AQ = 0, thus 



AYA"^H -f AQs = 0 



The node equations are then formulated 

YnH = Q 
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where 



Yn = AYA"^ 

and 

Q = -AQs 

Therefore, the node pressures H, can now be determine by 

h = y;'q 

Although many techniques exist for the solution of large sets of linear si- 
multaneous algebraic equations, the most efficient computationally, appears to be 
the Cholesky reduction followed by back subsitution that is employed in the Gauss 
elimination method. The only restriction on the use of Cholesky reduction is that 
its use is confined to symmetric, positive definite matrices. Here it is fortunate that 
Yn is always positive definite. 

The Cholesky reduction is based on the premise that there is a unique lower 
triangular matrix that permits a factorization in the form of A = if the matrix 
A is symmetric and positive definite. Computationally, the Cholesky reduction is a 
very efficient technique in that it only requires n(n -b l)/2 storage locations for L, 
rather than the usual v} locations required by other methods. [Ref. 3]. The number 
of operations using the Cholesky reduction is approximately n^/6 rather than the 
usual n^/3 required for most other decompositions [Ref. 3]. 

The Cholesky method to solve the basic system 



AX = B 



is based on finding the solution to two equivalent systems: 

L'^C = B and LX = C 
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The elements of C are determined by the algorithms 







and 



Ct = 



1-1 

£=1 

Sii 



i > 1 



Once C is known, X can be found using back substitution as employed in the Gauss 



elimination method, [Ref. 5] that is 



X 



n 



Cn 



S 



nn 



and 



C| 



X, = 



n 

^t£^l 

£=t-^l 



i < n 
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III. DESCRIPTION OF THE COMPUTER 

CODE 



A. PROGRAM STRUCTURES AND CAPABILITIES 

1. Restrictions and Limitations 

The computer system used is an IBM or IBM compatible personal com- 
puter. The current program fixes the maximum number of individual fiuid passages 
allowed at 300. For proper use of this program, the computer must have a minimum 
storage requirement of 2 M Bytes. A detailed computing time study for the program 
has not been undertaken because the computing time changes exponentially with 
the number of nodes and branches selected in the proposed cooling rack design. 

The computer code can be operated using the following information 

• Laminar Flow Conditions 

• Density of water varying from x to y 

• Viscosity of water varying from x to y 

• Up to 100 branches 

• Up to 100 nodes 

• Up to 40 sources 

• Any diameter (or equivalent diameter) flow passages 

• Branches of any length 

• Metric or English units 
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2. Structure of the Main Computer Code 

The the main computer code including the three associated subroutines 
are in Appendix A. The main code is essentially divided into two major sections, 
initialization/input and laminar flow solution with its associated output. The main 
function of the initialization/input section are: 

1 . Set up the problem formulation by reading in the necessary values (node from, 
node to, length and cross sectional dimensions) for each branch of the cooling 
rack. 

2. Read in the viscosity and density values. 

3. Read in location and characteristics of each pressure source. 

The second section calculates the pressure and flow distribution in the 
rack. It includes the subroutines MATMUL, DECOMP and CHOLESKY. The main 
functions of this section are to 

1. Develop the node-branch incidence matrix. A, using the information read in 
section one. 

2. Calculate the effective branch lengths to account for losses at the nodes. 

3. Calculate the area and effective diameter of each flow passage cross sectional 
area. 

4. Calculate the branch admittance matrix, Y. 

5. Develop the flow source matrix, Qg, from information read in section one. 

6. Calculate the transpose of the matrix A. 
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7. Calculate the node admittance matrix, Yq. 

8. Calculate the node flow source vector, Is- 

9. Calculate the node pressure vector, H. 

10. Calculate the inverse of Yn- 

11. Calculate the branch pressure vector, P. 

12. Calculate the branch flow rate vector, Q. 

13. Calculate the branch Reynolds numbers. 

14. Provide readouts of all branch flow rates, all node pressures and all branch 
Reynolds numbers. 

B. DESCRIPTION OF SUBROUTINES 

1. Subroutine MATMUL 

This subroutine [Ref. 5] multiplies an n x m matrix by any mxi matrix 
to form an n X / matrix and is called whenever matrix multiplication is required. 

2. Subroutine DECOMP 

This subroutine [Ref. 6] performs the decomposition of the symmetric, 
positive definite matrix, Y,i using the Cholesky reduction. It is called by subroutine 
CHOLESKY. 

3. Subroutine CHOLESKY 

This subroutine [Ref. 7] provides the solution of a linear system of 
equations using the Cholesky decomposition method for positive definite matrices. 
This subroutine calls subroutine DECOMP. 
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IV. RESULTS AND DISCUSSION OF CASE 

RUNS 

A. CASE RUN ONE 

Case one was run using the configuration shown in Figure 2.1 and the follow- 
ing data; 

1. Water density = 62.4 Ib/ft^ 

2. Water viscosity = 8.8 x 10““* Ib/sec-ft 

3. Basic shoebox design with 8 nodes and 12 branches 

4. Length of end branches (1-4 and 9-12) = 2.2 ft 

5. Length of middle branches (5-8) = 3.0 ft 

6. Circular passages with a diameter of 0.3 inches 

7. One flow source at branch one 

8. Strength of source was 0.4 gal/min 

The results from case one are in Appendix B. It should be noted that the 
Reynolds number in each branch indicates laminar flow and that there is cooling 
water in all sections of the sample rack design. These results demonstrate that 
the analysis technique presented here, does allow for the determination of whether 
or not a proposed cooling rack design will indeed provide a proper distribution of 
cooling water. It also verifies that the laminar flow assumption is valid for some 
design parameters. 
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B. CASE RUN TWO 

Case run two demonstrates the flexibility of the program through a varia- 
tion in some of the input parameters. The basic configuration is still the same 
configuration shown in Figure 2.1. 

1. Water density = 999.5 kg/m^ 

2. Water viscosity = 13.1 x 10“^ kg/m-sec 

3. Basic shoebox design with 8 nodes and 12 branches 

4. Length of end branches (1-4 and 9-12) = 1.0 m 

5. Length of middle branches (5-8) = 1.2 m 

6. Rectangular passages with a height of 1.3 cm and width of 1.4 cm. 

7. One flow source at branch one 

8. Strength of source was 1.8 lit/min 

The results from case two are located in Appendix B. 
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V. CONCLUSION 



A. GENERAL COMMENTS 

This computer code has been developed as an initial step in the total design 
of a modular cooling rack for avionics equipment. In itself, the code details a 
specific design technique and allows for the determination of whether one proposed 
configuration, including source location, characteristics of the cooling water, and the 
size and shape of the proposed flow passages will indeed provide a proper distribution 
of the cooling water. Without proper coolant distribution, the cooling rack will be 
inefficient and perhaps, totally ineffective. 

B. ENHANCING THE MAIN PROGRAM 
CAPABILITY 

The next step in the development of a computer code to provide more flex- 
ibility and range in the rack design is the inclusion of the turbulent flow regime 
within the coolant passages. After turbulent flow is effectively incorporated, the 
final phase in the total rack design is a modification for heat input in the individual 
coolant flow passages. 
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APPENDIX A 



PROGRAM COOL 

c Cooling Rack Calculations Using Node Analysis 

c INITIALIZATION 

r MU: fluid viscosity (Ib/ft-sec) 

:: RHO: fluid density {lb/ft"'3) 

r N: number of junctions in system 

:: N1 ; number of junction minus the datum node 

r B: number of branches in system 

:: G; acceleration of gravity (ft/sec'^2) 

r PI : pi 

INTEGER N,N1,B,I, J,K,QB,S 
INTEGER NT(40) ,NF(40) ,X, lUNITS 

REAL MU,RHO,G,PI,A(40,40) ,D(40) ,ELL(40) ,R(40) , MUl 
REAL QS (40,1), AT (40, 40) , QSl , I 1 (4 0 , 1 ) , LENGTH 
REAL E(40,1),P(40,1) ,Y(40,40) , YNI (4 0,40) 

REAL YP(40,1),Q(40,1),AY(40,40),YN(40,40) 

REAL V(40,l),RE(40,l),al(40),bl(40) , AREA (4 0) 

REAL IS (40, 1) , ISI (40, 1) ,RH01 
CHARACTER ANS , ANSI , ANS2 



G=32 . 174 
PI=3 . 1415926 
IUNITS=0 






THIS PROGRAM DETERMINES THE FLOW RATE AND PRESSURE DISTRIBUTION 
OF COOLING WATER WITHIN A VARIABLE -SI ZED AVIONICS COOLING RACK. 
VARIABLES INCLUDE RACK DIMENSIONS AND WATER CHARACTERISTICS. 



■ 2 * ******* -k * * -k ************************************************ -k * -k ***** -k * * * * 



1040 

1042 

1043 



FORMAT (///' ENTER THE TOTAL NUMBER OF NODES IN THE RACK SYSTEM 
' ,2X,\) 

FORMAT (BN, 13) 

FORMAT (//' ENTER THE TOTAL NUMBER OF BRANCHES IN THE RACK SYSTEM 
' ,2X,\) 



222 WRITE (lOT, 1040) 

READ (IN, 1042) N 
WRITE (lOT, 1045) N 

IF(N .LT. 0 .OR. N . GT . 100) THEN 
WRITE (lOT, 6565) 

GO TO 222 
END IF 

223 WRITE (lOT, 1043) 

READ (IN, 1042) B 
WRITEdOT, 1046) B 

IF(B .LT. 0 .OR. B . GT . 100) THEN 
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WRITE (lOT, 6567) 
GO TO 223 
END IF 



045 


FORMAT 


(//' 


YOU ENTERED ' 


046 


FORMAT 


(//' 


YOU ENTERED ' 


565 


FORMAT 


(/' 


THE NUMBER OF 


+ 


IX, THAN 


100 . ' 


A) 


567 


FORMAT 


(/' 


THE NUMBER OF 



THE NUMBER OF BRANCHES MUST BE GREATER THAN ZERO AND LESS 



+ THAN 100 . ' , 



N1=N-1 
DO 5 1=1, B 
QS(I,1)=0. 
5 CONTINUE 



5 CONTINUE 

WRITE (lOT, 7676) 

576 FORMAT (/' Do you want to work in British units or SI units (B/S)? 

+ ' A) 

READ (IN, 7677) ANS 

577 FORMAT (Al) 

IF(ANS .EQ. 'B' .OR. ANS .EQ. 'b') GO TO 7800 
IF(ANS .EQ. 'S' .OR. ANS .EQ. 'S') GO TO 7678 
GO TO 55 

300 lUNITS = 1 

WRITE (lOT, 1000) 

DOO FORMAT (//' Input Viscosity ( x 10"^4 Ibm/f t -sec) ' , 2X, \) 

READ (IN, 1003) MU 
MU1=MU* . 0001 
WRITE (lOT, 1002) 

D02 FORMAT (/' Input density (Ib/f t^3 ) ' , 2X, \) 

READ (IN, 1001) RHO 

301 FORMAT (BN, F8.4) 

303 FORMAT (BN, F7.6) 

304 FORMAT (BN, 13) 

GO TO 7801 



578 WRITE (lOT, 7602) 

502 FORMAT (//' Input Viscosity (Kg/m-sec x 10^4)',2X,\) 
READ (IN, 1003) MU 

WRITE (lOT, 7603) 

503 FORMAT (/' Input Density (Kg/m^3 ) ' , 2X, \) 

READ (IN, 1001) RHOl 

MU1=MU* . 00006719 
RHO=RH01* .06243 



A: NODE-BRANCH INCIDENCE MATRIX 
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c INITIALIZE A MATRIX 
c 

7801 DO 20 1=1, N1 

DO 20 J=1,B 

20 A(I,J)=0. 

c 

DO 22 1=1, B 
NF (I) =0 
22 NT(I)=0 

c 
c 

c DETERMININATION OF PASSAGE SHAPE 
c 
c 
c 

8000 WRITE (lOT, 8001) 

8001 FORMAT (/' IF FLOW PASSAGES ARE CIRCULAR, ENTER A ZERO.',/, 

+ ' IF PASSAGES ARE RECTANGULAR, ENTER A ONE.',2X,\) 

READdN, 1042) X 

IF(X .EQ. 0) THEN 
GO TO 8004 
END IF 

IF(X .EQ. 1) THEN 
GO TO 8005 
END IF 

8003 WRITE (lOT, 8006) 

8006 FORMAT (/' ERRONEOUS INPUT' \) 

GO TO 8000 
c 
c 

c INITIAL DATA INPUT TO DEVELOP A, ELL AND D MATRICES FOR 
c RECTANGULAR PASSAGES 
c 

c ELL : LENGTH 

c D: DIAMETER 

c 
c 

8005 CALL CLS 

IF (lUNITS .EQ. 0) THEN 
DO 1493 J=1,B 
WRITE (lOT, 1051) J 
WRITE (lOT, 1492) 

1492 FORMAT(/' FROM NODE, TO NODE, LENGTH (m) , HEIGHT(cm), WIDTH(cm)', 

+ /) 

READdN, 8202) NF(J) ,NT(J) ,ELL(J) ,al(J) ,bl(J) 

1493 CONTINUE 

7272 FORMAT (/' YOU ENTERED THE FOLLOWING DATA.'\) 

7273 FORMAT (/' BRANCH FROM NODE TO NODE LENGTH (m) HEIGHT (cm) ', 2X, 

+ 'WIDTH (cm) ' , /) 

7274 FORMAT (2X, 13 , 6X, 13 , 8X, 13 , 2X , F9 . 5 , 4X , F7 . 5 , 4X , F7 . 5 ) 

7275 FORMAT (/' IS THIS CORRECT (Y/N) ? ' \) 

CALL CLS 

16 WRITE (lOT, 7272) 

WRITE (lOT, 7273) 

DO 7277 J=1,B 

WRITE (lOT, 7274) J, NF ( J) , NT ( J) ,ELL(J) ,al(J) ,bl(J) 
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CONTINUE 
WRITE (lOT, 7275) 

READ (IN, 7677) ANS 

IF(ANS .EQ. 'Y' .OR. ANS .EQ.'y') GO TO 7278 

IF(ANS .EQ. 'N' .OR. ANS . EQ . 'n') GO TO 7280 

GO TO 7276 

DO 7279 J=1,B 

ELL(J) =ELL(J) *3.28 

al (J) =al (J) * .3937 

bl(J)=bl(J)*.3937 

CONTINUE 

GO TO 1495 

CONTINUE 

FORMAT ( / ' INPUT ONE BRANCH NUMBER TO CHANGE VALUES . ' \ ) 
WRITE (lOT, 7281) 

READ (IN, 1042) K1 

FORMAT (/' INPUT VALUES FOR BRANCH IX, 13 , \) 

WRITE (lOT, 7282) K1 
WRITE (lOT, 1492) 

READ (IN, 8202)NF(K1) ,NT(K1) ,ELL(K1) , al (K1 ) , bl (K1 ) 

FORMAT (/' DO YOU HAVE ANYMORE CHANGES (Y/N)?'\) 

WRITE (lOT, 7283) 

READ (IN, 7677) ANSI 

IF (ANSI .EQ. 'Y' .OR. ANSI . EQ . 'y') GO TO 16 
IF (ANSI .EQ. 'N' .OR. ANSI . EQ . 'n') GO TO 7278 
GO TO 7285 
END IF 

CONTINUE 



DO 8200 1=1, B 
WRITE (lOT, 1051) I 
WRITE (lOT, 8201) 

FORMAT(/' FROM NODE, TO NODE, LENGTH(FT), HEIGHT(IN), WIDTH(IN) 

' ,/) 

FORMAT (2 13 , F7 . 4 , 2F6 . 4 ) 

READ (IN, 8202) NF(I) ,NT(I) ,ELL(I) ,al(I) ,bl(I) 

CONTINUE 

CALL CLS 
WRITE (lOT, 7272) 

WRITE (lOT, 8207) 

FORMAT(/' BRANCH FROM NODE TO NODE LENGTH(FT) HEIGHT ( IN) ', 2X, 
'WIDTH (IN) ' ,/) 

FORMAT (2X, 13 , 7X, 13 , 7X, 13 , 4X, F9 . 5 , 6X , F7 . 5 , 3X , F7 . 5 ) 

DO 8277 J=1,B 

WRITE (lOT, 7374) J,NF(J) ,NT(J) ,ELL(J) , al (J) ,bl (J) 

CONTINUE 
WRITE (lOT, 7275) 

READ (IN, 7677) ANS 

IF(ANS .EQ. 'Y' .OR. ANS .EQ.'y') GO TO 1495 

IF(ANS .EQ. 'N' .OR. ANS . EQ . 'n') GO TO 8280 

GO TO 8276 

CONTINUE 

WRITE (lOT, 7281) 
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READ (IN, 1042) K1 
WRITE (lOT, 7282) K1 
WRITE (lOT, 8201) 

READ (IN, 8202)NF(K1) ,NT(K1) ,ELL(K1) , al (K1 ) , bl (Kl) 

8285 WRITE (lOT, 7283) 

READ (IN, 7677) ANSI 

IF (ANSI .EQ. 'Y' .OR. ANSI .EQ. 'y') GO TO 15 
IF (ANSI .EQ. 'N' .OR. ANSI .EQ. 'n') GO TO 1495 
GO TO 8285 

1495 DO 1496 1=1, B 

IF (NF(I) .LT. N) THEN 
A(NF(I) , I) = 1. 

END IF 

IF (NT (I) .LT. N) THEN 

A(NT(I) , I) = -1. 

END IF 
al (I) =al (I) /12 
bl (I) =bl (I) /12 
AREA(I) =al (I) *bl (I) 

D(I) = (2*(al(I)*bl(I) ) )/(al(I)+bl(I)) 

LENGTH=ELL(I) 

ADDL=60*d(I) 

ELL ( I ) = LENGTH + ADDL 

1496 CONTINUE 
GO TO 8300 

c 

c 

c RECEIVE DATA FROM KEYBOARD TO DEVELOP A,D AND ELL MATRICES FOR 
c CIRCULAR PASSAGES 
c 

1050 FORMAT (IX, 13 , 8X, 13 , 5X, F7 . 4 , 4X, F6 . 4 ) 

9050 FORMAT (213, 2F7 .4) 

8004 CALL CLS 

IFdUNITS .EQ. 0) THEN 

GO TO 993 

ELSE 

DO 25 1=1, B 

1051 FORMAT (/' INPUT THE FOLLOWING DATA, SEPARATED BY COMMAS FOR, IX, 
+ BRANCH: ' , 13, 2X, \) 

1053 FORMAT (/' FROM NODE , TO NODE , LENGTH (FT) , DIAMETER ( IN) ', 2X ,/ ) 

WRITE (lOT, 1051) I 
WRITE (lOT, 1053) 

READ (IN, 9050) NF(I) ,NT(I) ,ELL(I) ,D(I) 

25 CONTINUE 

CALL CLS 

18 WRITE (lOT, 7272) 

WRITE (lOT, 1059) 

1059 FORMAT (/' BRANCH FROM NODE TO NODE LENGTH (FT) DIAMETER ( IN) ' , 

+ / ) 

DO 9277 J=1,B 

WRITE (lOT, 1155) J,NF(J) ,NT(J) ,ELL(J) ,D(J) 

1155 FORMAT (2X, 13 , 9X, 13 , 8X, 13 , 6X, F7 . 4 , 4X, F6 . 4) 

1255 FORMAT (2X , 13 , 7X , 13 , 7X, 13 , 5X , F7 . 4 , 5X, F6 . 4 ) 

9277 CONTINUE 
9276 WRITE (lOT, 7275) 

READ (IN, 7677) ANS 
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IF(ANS .EQ. 'Y' .OR. ANS .EQ.'y') GO TO 26 
IF(ANS .EQ. 'N' .OR. ANS . EQ . 'n') GO TO 9280 

GO TO 9276 
280 CONTINUE 

284 WRITE (lOT, 7281) 

READdN, 1042) K1 
WRITE (lOT, 7282) K1 
WRITE (lOT, 1053) 

READdN, 9050)NF(K1) ,NT(K1) ,ELL(K1) ,D(K1) 

285 WRITE (lOT, 7283) 

READ (IN, 7677) ANSI 

IF(ANS1 .EQ. 'Y' .OR. ANSI .EQ. 'y') GO TO 18 
IF (ANSI .EQ. 'N' .OR. ANSI .EQ. 'n') GO TO 26 
GO TO 9285 
END IF 
GO TO 26 



?3 DO 29 J=1,B 

WRITE (lOT, 1051) J 

355 FORMAT (/' FROM NODE, TO NODE, LENGTH (M) , DIAMETER (CM) ',/ ) 

WRITE (lOT, 1055) 

READdN, 9050) NF(J) ,NT(J) ,ELL(J) ,D(J) 

3 CONTINUE 

CALL CLS 

3 WRITE (lOT, 7272) 

WRITE (lOT, 1077) 

377 FORMAT (/' BRANCH FROM NODE TO NODE LENGTH (M) DIAMETER (CM) ' , 

+ /) 

DO 1277 J=1,B 

WRITE (lOT, 1255) J,NF(J) ,NT(J) ,ELL(J) ,D(J) 

111 CONTINUE 

>76 WRITE (lOT, 7275) 

READ (IN, 7677) ANS 

IF (ANS .EQ. 'Y' .OR. ANS .EQ.'y') GO TO 111 
IF(ANS .EQ. 'N' .OR. ANS . EQ . 'n') GO TO 1280 

GO TO 1276 
>80 CONTINUE 

>84 WRITE (lOT, 7281) 

READdN, 1042) K1 
WRITE (lOT, 7282) K1 
WRITE (lOT, 1055) 

READdN, 9050)NF(K1) ,NT(K1) ,ELL(K1) ,D(K1) 
i:85 WRITE (lOT, 7283) 

READ (IN, 7677) ANSI 

IF (ANSI .EQ. 'Y' .OR. ANSI .EQ. 'y') GO TO 19 
IF (ANSI .EQ. 'N' .OR. ANSI .EQ. 'n') GO TO 111 
GO TO 1285 

II DO 28 J=1,B 

ELL(J) =ELL(J) *3.28 
D (J) =D (J) * .3937 
: CONTINUE 



DO 33 1=1, B 
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IF (NF(I) .LT. N) THEN 
A(NF(I) , I) =1. 

END IF 

IF (NT (I) .LT. N) THEN 
A(NT(I) ,!)=-!. 

END IF 

D(I) =D(I) /12 

AREA (I) =PI* (D (I) **2) /4 

LENGTH=ELL (I) 

ADDL=60*D (I) 

ELL ( I ) =LENGTH+ADDL 

33 CONTINUE 

c 
c 
c 
c 

8300 DO 30 K=1,B 

R(K) = (RHO*G*AREA(K) * ( (D(K) ) **2) ) / (32 . *MU1*ELL(K) ) 
30 CONTINUE 

c 

DO 35 J=1,B 

DO 35 K=1,B 
Y(J,K) =0. 

35 CONTINUE 

c 
c 

DO 40 1=1, B 

Y ( I , I ) =R ( I ) 

40 CONTINUE 



c 

c 



c QS: 


MATRIX 


OF FLOW SOURCES 








c 

C S: 


NUMBER 


OF FLOW SOURCES 








c QB: 


SOURCE 


BRANCH 








C QSl : 


SOURCE 


STRENGTH 








c 

1015 


FORMAT 


(/' Input the number of 


flow sources;' 


,3X,\) 


1016 


FORMAT 


(//' Input the branch of 


source ' , 13 , 3X 


, A) 


1017 


FORMAT 


(//' Input the strength 


of 


the source 


(gal/min) : ' , 3X, \) 


9017 


FORMAT 


(//' Input the strength 


of 


the source 


(liter/min) : ' , 3X, /) 



WRITE (lOT, 1015) 

READdN, 1004) S 
DO 50 1=1, S 

WRITE (lOT, 1016) I 
READdN, 1004) QB 

IFdUNITS .EQ. 1) THEN 
WRITE (lOT, 1017) 

READ (IN, 1001) QSl 
ELSE 

WRITE (lOT, 9017) 
READdN, 1001) QSl 
QS1=QS1* .2642 
END IF 
QSl=QSl/444 
DO 60 K=1,B 
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IF (K .EQ. QB) THEN 

QS (K, 1) =QS1 

ELSE 

QS(K,1)=0. 

END IF 

50 CONTINUE 

50 CONTINUE 



WRITE (lOT, 3012) (QS ( I , 1) , 1=1 , B) 



>1ATRIX MANIPULATION TO DETERMINE INDIVIDUAL BRT^CH FLOW AND 
PRESSURE DISTRIBUTIONS 

AT: TRANSPOSE OF MATRIX A 

AY: PRODUCT OF MATRIX A AND MATRIX Y 

YN: PRODUCT OF MATRIX AY AND MATRIX AT (YN=AYAT) 



DO 80 1=1, N1 
DO 80 J=1,B 
30 AT(J, I) =A(I, J) 



CALL MATMUL(AY,A,Y,N1,B,B) 
CALL MATMUL(YN,AY,AT,N1,B,N1) 



IS: NODE FLOW SOURCE VECTOR IS= -AQS 
CALL MATMUL(IS,A,QS,N1,B, 1) 



DO 90 1=1, N1 

IS(I,1)=-IS(I,1) 
ISI(I,1)=IS(I,1) 
10 CONTINUE 



YNI: MATRIX EQUAL TO MATRIX YN BEFORE CHOLESKI INVERSION. AFTER 
DECOMPOSITION IT HOLDS THE UPPER TRIANGULAR MATRIX. 

ISI: MATRIX EQUAL TO MATRIX IS BEFORE INVERSION. AFTER INVERSION, 
IT HOLDS THE SOLUTION TO E=YN(-l)IS 



DO 95 1=1, N1 
DO 95 J=1,N1 

YNI (I, J)=YN(I, J) 
CONTINUE 
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CALL CHOLESKY(YNI, ISI,N1,N1) 
c 
c 

C E: NODE VOLTAGE MATRIX E=YN(-1)IS 
c P: BRANCH PRESSURE P=ATE 

c Q: BRANCH FLOW RATE Q=YP+QS 

c YP: PRODUCT OF MATRIX Y AND MATRIX P (YP) 
c 

DO 96 1=1, N1 

E (I, 1) =ISI (1,1) 

96 CONTINUE 

c 
c 

CALL MATMUL(P,AT,E,B,N1,1) 
c 
c 

CALL MATMUL(YP,Y,P,B,B,1) 
c 

DO 120 1=1, B 

Q(I,1)=YP(I,1)+QS(I,1) 

120 CONTINUE 

CALL CLS 
c 
c 

c REYNOLDS NUMBER CALCULATIONS 
c 

DO 5000 1=1, B 

V{I,1)=ABS(Q(I,1)/(AREA{I) ) ) 

RE ( I , 1 ) =rho* V ( I , 1 ) *d { I ) /mul 
IF(RE(I,1) .GT. 2100) THEN 
WRITE (lOT, 5001) I 
WRITE (lOT, 5002) 

5001 FORMAT(/' REYNOLDS NUMBER IN BRANCH' , 13 , 2x, ' EXCEEDS ', \) 

5002 FORMAT (/' 2100. THE FLOW IS NOT LAMINAR 2X, \) 

ELSE 

END IF 

5000 CONTINUE 

c 
c 

1500 FORMAT (3X, 13 , 7X, Fll . 8 , 2X, 2F15 . 5 , IX, FIO .0) 

1700 FORMAT (3X, 13 , 5X, Fll . 8 , 3X, 2F15 . 5 , 2X, F5 . 0) 

1501 FORMAT{///' Branch' , 4X, ' Flow Rate (ft/sec) ', 4x, ' P in (psf)',6x, 

+ 'P out (psf ) ' , 6x, ' Re' , lx, //) 

1601 FORMAT (///' Branch' , 4x, ' Flow Rate (m/sec) ', 5x, ' P in (N/m^2)',7x, 
+ 'P out (N/m^2) ' , 6x, ' Re' , 2x, //) 

IFdUNITS .EQ. 1) THEN 
WRITE (lOT, 1501) 

DO 7000 1=1, B 

WRITE(IOT,1500) I,Q(I,1) , E ( (NF ( I ) ) , 1 ) , E ( (NT ( I ) ) , 1 ) , RE ( I , 1 ) 

7000 CONTINUE 

99 FORMAT(/' WOULD YOU LIKE A PRINTOUT OF THIS TABLE (Y/N)?'\) 

98 WRITE (lOT, 99) 

READ (IN, 7677) ANS2 

IF(ANS2 .EQ. 'Y' .OR. ANS2 . EQ . 'y') GO TO 97 
IF(ANS2 .EQ. 'N' .OR. ANS2 . EQ . 'n') GO TO 83 
GO TO 98 
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WRITEdPR, 1501) 

DO 93 1=1, B 

WRITEdPR, 1500) 1,0(1,1) ,E( (NF(I) ) ,1) ,E( (NT(I) ) ,1) ,RE(I,1) 

CONTINUE 

CONTINUE 

ELSE 

DO 7470 J=1,N1 
E(J,1)=E(J,1)*47.88 
70 CONTINUE 

WRITE (lOT, 1601) 

DO 7001 1=1, B 
Qd,l)=Q(I,l)*.3048 

WRITE (lOT, 1500) I , Q (I , 1) , E ( (NF (I ) ) , 1) , E ( (NT (I ) ) , 1) , RE (I , 1) 
01 CONTINUE 

WRITE (lOT, 99) 

READ (IN, 7677) ANS2 

IF(ANS2 .EQ. 'Y' .OR. ANS2 .EQ. 'y') GO TO 47 
IF(ANS2 .EQ. 'N' .OR. ANS2 .EQ. 'n') GO TO 82 
GO TO 48 
WRITE (IPR, 1601) 

DO 43 1=1, B 

WRITE (IPR, 1500) I, Q (1,1) , E ( (NF (I ) ) , 1 ) , E ( (NT (I ) ) , 1 ) , RE (I , 1 ) 

CONTINUE 

CONTINUE 

END IF 

END 
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SUBROUTINE CHOLESKY (A, B , N, NX) 



c 

c 

c SOLUTION OF LINEAR SYSTEMS OF EQUATIONS BY THE CHOLESKI 
c METHOD FOR SYMMETRIC POSITIVE DEFINITIVE MATRICES, 
c 

c A: ARRAY CONTAINING THE SYSTEM MATRIX (AX=B) 

c B: ARRAY CONTAINING THE VECTOR IF INDEPENDENT COEFFICIENTS 

c C: AUXILIARY VECTOR 

c N: ORDER OF A 

c NX: ROW AND COLUMN DIMENSION OF A 
c 

DOUBLE PRECISION B(NX,1) 

DIMENSION A(40,40) 
c 

c COMPUTE UPPER TRIANGULAR MATRIX FROM A AND STORE ALSO IN A 



CALL DECOMP (A, N, NX) 
c 

c COMPUTE THE C VECTOR AND STORE IN ARRAY B 
c 

B(1,1)=B(1,1)/A(1,1) 





WRITE (lOT, 100) 


A(l, 1) 


100 


FORMAT (F18. 15) 






DO 10 1=2, N 






WRITE (lOT, 101) 


I,N 


101 


FORMAT (213) 






WRITE (lOT, 102) 


A(I, I) ,A(N,N) 


102 


FORMAT (2F13 .11) 





D=B(I, 1) 
11 = 1-1 



DO 5 L=l, II 

5 D=D-A(L, I) *B (L, 1) 

10 B(I,1)=D/A(I,I) 

B(N,1)=B(N,1)/A(N,N) 

c 

c COMPUTE THE SYSTEM UNKNOWNS AND STORE IN ARRAY B 



c 

N1=N-1 

DO 30 L=1,N1 
K=N-L 

C WRITE (lOT, 105) A(K,K),K 

105 FORMAT (F13.il, 13) 

K1=K+1 

DO 20 J=K1,N 

20 B (K, 1) =B (K, 1) -A(K, J) *B (J, 1) 

30 B(K,1)=B(K,1)/A(K,K) 

RETURN 

END 



36 



SUBROUTINE MATMUL (C, A, B, N, M, L) 

This subroutine computes the matrix operation C = A * B 

N: NUMBER OF ROWS IN MATRIX A AND C 
M: NUMBER OF COLUMNS IN A AND ROWS IN B 
L: NUMBER OF COLUMNS IN B AND C 

DIMENSION A(40,40),B(40,40),C(40,40) 

DO 20 1=1, N 
DO 20 J=1,L 
C(I, J) =0.0 
DO 20 K=1,M 

) C(I, J)=C(I,J)+A(I,K) *B(K, J) 

RETURN 

END 



37 



c 



SUBROUTINE DECOMP (A, N, NX) 



c THIS SUBROUTINE PERFORMS THE DECOMPOSITION OF A POSITIVE, 
c DEFINITE, SYMMETRIC MATRIX INTO AN UPPER TRIANGULAR MATRIX, 
c 

c A: ARRAY ORIGINALLY CONTAINING MATRIX TO BE DECOMPOSED, 
c AT THE END IT CONTAINS THE UPPER TRIANGULAR MATRIX, 

c N: ORDER OF A 

c NX: ROW AND COLUMN DIMENSION OF A 

INTEGER N 
DIMENSION A(40,40) 

IF (A(l,l) )1,1,3 

1 WRITE(IOT,2) 

2 FORMAT (' ZERO OR NEGATIVE RADICAND' ) 

C WRITE (lOT, 500) A(l,l) 

500 FORMAT (' A ( 1 , 1 ) ' , F13 . 11 ) 

GO TO 200 

3 A(1,1)=SQRT(A(1,1) ) 

DO 10 J=2,N 

10 A(l, J) =A(1, J) /A(l, 1) 

DO 40 I=2,N 

C WRITE (lOT, 600) I 

600 FORMAT (13) 

11 = 1-1 
D=A(I, I) 

DO 20 L=1,I1 

20 D=D-A(L, I) *A(L, I) 

C WRITE (lOT, 503) I,D,A(I,I) 

503 FORMAT ( 13, 2F13. 11) 

IF (A(I, I) ) 1, 1, 21 

21 A(I,I) =SQRT(D) 

IF (I .EQ. N) THEN 
GO TO 45 

END IF 

12 = 1+1 

DO 40 J=I2,N 
D=A(I, J) 

DO 30 L=1,I1 
30 D=D-A(L, I) *A(L, J) 

C WRITE (lOT, 503) I,D,A(I,I) 

4 0 A(I, J) =D/A(I, I) 

45 DO 50 1=2, N 

11 = 1-1 

DO 50 J=l, II 
50 A(I,J)=0 

200 RETURN 
END 
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2 

3 
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6 
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8 

9 

10 

11 

12 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 



APl’ICNDIX I) 



Case Run One 



Flow Rate 


Pin 


Pout 


Re 


(ft/sec) 


(psf ) 


(psf) 




.00035966 


- .03996 


. 05160 


1299 . 


. 00019242 


.05160 


. 01905 


695 . 


. 00015641 


. 01905 


- . 00741 


565 . 


. 00019242 


- .00741 


- . 03996 


695 . 


. 00016724 


- . 00555 


- . 03996 


604 . 


. 00016724 


. 05160 


.01719 


604 . 


. 00003601 


. 01905 


. 01164 


130 . 


. 00003601 


- . 00741 


. 00000 


130 . 


. 00013443 


- . 00555 


.01719 


485 . 


. 00003281 


. 01719 


. 01164 


118 . 


. 00006882 


. 01164 


. 00000 


249 . 


. 00003281 


. 00000 


- . 00555 


118 . 



Case Run Two 



Flow Rate 


Pin 


Pout 


Re 


(m/sec) 


(N/sm) 


(N/sm) 




. 00013296 


- .28579 


.37451 


698 


. 00006895 


.37451 


. 13925 


362 


. 00005561 


. 13925 


- . 05052 


292 


.00006895 


- . 05052 


- . 28579 


362 


. 00006401 


- .04323 


- . 28579 


336 


. 00006401 


.37451 


. 13195 


336 


. 00001333 


. 13925 


. 08872 


70 


.00001333 


- . 05052 


. 00000 


70 


.00005134 


- . 04323 


. 13195 


270 


.00001267 


. 13195 


. 08872 


67 


.00002600 


. 08872 


. 00000 


137 


.00001267 


. 00000 


- . 04323 


67 
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